###########################################
# Figure A1: SOI missing zip codes

###########################################

ca_zip <- zctas(starts_with = c("90","91","92","93","94","95","96"),year = 2016)
ca_zip <- ca_zip[as.numeric(substr(ca_zip$ZCTA5CE10,1,3))<=961,]
soi <- read.csv("Data/Raw data/soi_zipcode_panel_with_all_noagi.csv")

soi <- soi[soi$year == 2016,] # 2016 is the year with the most data

ca_zip$`Zip codes` <- ifelse(ca_zip$ZCTA5CE10 %in% soi$zip,"In SOI","Not in SOI")

ca <- states()
ca <- ca[ca$STUSPS == "CA",]

ca$`Zip codes` <- " "

p <- ggplot(data = ca_zip) + 
  geom_sf(data=ca,aes(fill = `Zip codes`))+
  geom_sf(aes(fill = `Zip codes`), linewidth = .1, color = "black")+
  scale_fill_manual(values = c(" " = "white", "Not in SOI" =  "#B31B1B","In SOI" = "#FFCCCC"), name = "Zip codes:")+
  scale_x_continuous()+
  theme_map()
p
ggsave(plot = p, file="Output/fig A1.Missing zipcodes in CA SOi.png" ,bg = "white")



rm(list = setdiff(ls(),c("df","city")))


###########################################
# Figure A2: Share of tax returns with property tax deduction

###########################################

soi <- read.csv("Data/Raw data/soi_school_districts.csv")
soi$School_District <- sprintf("%07d",soi$School_District)
soi <- soi[substr(soi$School_District,1,2)=="06",]#Keep California

# Aggregate by year the nb of itemizers and number of tax returns (mean, 5th and 95th quantiles)
b <- aggregate(x = list(itm_perc = soi$itm_perc),by = list(year = soi$year),FUN = mean)
b <- aggregate(x = list(itm_perc = soi$itm_prop_perc),by = list(year = soi$year),FUN = 'quantile', probs=c(0.05,0.5,0.95) )
b <- as.data.frame(cbind(b$year,b$itm_perc))
colnames(b) <- c("year","q05","q50","q95")
b$col <- "Color"
b$year  <- as.Date(paste(b$year,"06","30",sep = "-"))# set date in middle of the year for aesthetic

p <- ggplot(data = b,aes(x = year,y = q50,col = col))+
  geom_line()+
  geom_point()+
  scale_y_continuous(name = "Share of property tax deducters",labels = percent,limits = c(0,.6))+
  scale_x_date(name = "") +
  geom_ribbon(aes(ymin=q05, ymax=q95), linetype=3, alpha=0.1) +
  scale_color_manual(values=c("#B31B1B","#222222","#222222")) +
  geom_vline(xintercept= as.Date("2018-01-01"), linetype=5, color = "#222222") + 
  theme_minimal()+
  theme(legend.position = "none")

ggsave(plot = p, file="Output/fig A2.School_district_distribution_of_itmshare.png" ,bg = "white",width = 10,height = 5)

rm(list = setdiff(ls(),c("df","city")))

###########################################
# Figure A3: Intention to approve school district bonds and parcel levies

###########################################

df_ppic <- read.csv("Data/Raw data/ppic.csv")

df_ppic <- df_ppic[!is.na(df_ppic$year),]
# Removes NA and "do not know" responses
df_ppic <- df_ppic[!is.na(df_ppic$local_bond_acceptance),]
df_ppic <- df_ppic[!is.na(df_ppic$local_prop_or_parcel_tax_acceptance),]
df_ppic <- df_ppic[df_ppic$local_bond_acceptance != "dunno",]
df_ppic <- df_ppic[df_ppic$local_prop_or_parcel_tax_acceptance != "dunno",]

df_ppic$votes <- ifelse(df_ppic$local_bond_acceptance == "yes",1,0)
df_ppic$votes2 <- ifelse(df_ppic$local_prop_or_parcel_tax_acceptance == "yes",1,0)
df_ppic <- df_ppic[df_ppic$tenure != "dunno",]

# Keep voters 
df_ppic <- df_ppic[!is.na(df_ppic$voter),]

#######
# Obtain mean and standard errors by year

R <- lm(data = df_ppic,formula =  votes ~ 0 + as.factor(year),weights = df_ppic$weight)
a <- as.data.frame(cbind(R$coefficients,confint(R,level = 0.95)))
colnames(a) <- c("perc_yes","lb","ub") 
a$year <- as.numeric(substr(rownames(a),16,19))

R <- lm(data = df_ppic,formula =  votes2 ~ 0 + as.factor(year),weights = df_ppic$weight)
b <- as.data.frame(cbind(R$coefficients,confint(R,level = 0.95)))
colnames(b) <- c("perc_yes","lb","ub") 
b$year <- as.numeric(substr(rownames(b),16,19))

b$`Approval of` <- "Bonds"
a$`Approval of` <- "Parcel Levy"

all <- rbind(a,b)

#Plot
q <- ggplot(data = all,aes(x = year,y = perc_yes,group= `Approval of`,col = `Approval of`))+
  geom_line()+
  geom_point()+
  geom_ribbon(aes(ymin=lb, ymax=ub), alpha = 0.1,fill = "grey",colour = NA) +
  scale_y_continuous(name = "Intention to vote in favor of local public goods",labels = percent)+
  geom_vline(xintercept = 2018.5, linetype = "dashed")+
  scale_x_continuous(name = "")+
  scale_color_manual(values = c("black","#B31B1B"))+
  theme_minimal()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())


ggsave(plot = q, file="Output/fig A3.survey_approval_rates.png" ,bg = "white",width = 10,height = 5)

rm(list = setdiff(ls(),c("df","city")))


###########################################
# Figure A4: Frequency of Referendums by School Districts

###########################################

a <- as.data.frame(table(df$LEAID))
colnames(a) <- c("LEAID","Nb.of.ref")
df <- merge(df,a,by = "LEAID")

df$Nb.of.ref <- as.character(df$Nb.of.ref)
df$Nb.of.ref <- ifelse(df$Nb.of.ref %in% c(8,9,10),">=8",df$Nb.of.ref)

q <- ggplot(data = df,aes(x = Nb.of.ref,fill ="#B31B1B"))+
  geom_bar()+
  scale_fill_manual(values=c("#B31B1B")) + 
  scale_x_discrete(limits = c("1","2","3","4","5","6","7",">=8"))+
  labs(title = "",y = "Total number of referendums",x = "Number of referendums by school districts")+
  theme_minimal()+
  theme(legend.position = "none")


ggsave(plot = q, file="Output/fig A4.Frequency_of_referenda.png" ,bg = "white",width = 10,height = 5)

rm(list = setdiff(ls(),c("df","city")))


###########################################
# Figure A5: Change in ``Yes'' votes and treatment intensity

###########################################

a <- aggregate(df$perc_yes,by = list(LEAID = df$LEAID,post = df$post),FUN = mean)
a_pre <- a[a$post==0,]
a_pre$post <- NULL
a_post <- a[a$post ==1,]
a_post$post <- NULL
colnames(a_pre)  <- c("LEAID","mean_pre")
colnames(a_post)  <- c("LEAID","mean_post")
a <- merge(a_pre,a_post,by = "LEAID")

a$diff_outcome <- a$mean_post-a$mean_pre

b <- df[,c("LEAID","itm_prop_change")]
a <- merge(a,b,by="LEAID",all.x = TRUE)

# Define groups 
groups_cut <- c(0.0,0.075,.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.3)

a2 <- aggregate(a, #the data frame
                by=list(cut(a$itm_prop_change,groups_cut)), #the bins (see below)
                FUN = mean) #the aggregating function

a3 <- aggregate(a$LEAID, #the data frame
                by=list(cut(a$itm_prop_change,groups_cut)), #the bins (see below)
                FUN = length) #the aggregating function

a2 <- merge(a2,a3,by = "Group.1")


q <- ggplot(a2, aes(x=itm_prop_change*100, y=diff_outcome, size=x)) +
  geom_point(alpha=0.75,color = "#B31B1B") +
  scale_size(range=c(1, 15))+
  geom_hline(yintercept=0, linetype=5, color = "#91959C") + 
  scale_x_continuous(name = "Chg.Ded (pp)",labels =  number,limits = c(0,25))+
  scale_y_continuous(name= "Difference in Yes votes (pp)",labels =  number,limits = c(-12,0))+
  theme_minimal()+
  theme(legend.position="none")

ggsave(plot = q, file="Output/fig A5.Empirical_ATT.png" ,bg = "white",width = 8,height = 5)
rm(list = setdiff(ls(),c("df","city")))


###########################################
# Figure A6: Distribution of the Treatment Intensity Variable

###########################################

###########################################

q <- ggplot(data = df,aes(x = itm_prop_change))+
  geom_histogram(aes(y = ..density..),
                 colour = 1, fill = "white",bins = 30)+
  geom_density(color = "#B31B1B",linewidth =1.5)+
  scale_x_continuous(name = "Change.Ded",labels = percent)+
  theme_minimal()

ggsave(plot = q, file="Output/fig A6.Ded_density.png" ,bg = "white")
rm(list = setdiff(ls(),c("df","city")))


###########################################
# Figure A7: California Local Governments Revenue Sources

###########################################

# All definitions of the Annual Survey of State and Local Government Finances can be found here: https://www.census.gov/programs-surveys/gov-finances.html

slgf <- read.csv("Data/Raw data/SLGF_2017.csv",header = TRUE)

# Keep California
slgf <- slgf[slgf$state ==5,]


list_exp <- c("E01","E03","E04","E05","E12","E16","E18","E21","E22","E23","E24","E25","E26","E29","E31","E32","E36","E44","E44","E45","E50","E52","E55","E56","E59","E60","E61","E62","E66","E74","E75","E77","E79","E80","E81","E85","E87","E89","E90","E91","E92","E93","E94","I89","I91","I92","I93","I94","J19","L67","J68","J85","X11","X12","Y05","Y06","Y14","Y53","F01","F03","F04","F05","F12","F16","F18","F21","F22","F23","F24","F25","F26","F29","F31","F32","F36","F44","F45","F50","F52","F55","F56","F59","F60","F61","F62","F66","F77","F79","F80","F81","F85","F87","F89","F90","F91","F92","F93","F94","G01","G03","G04","G05","G12","G16","G18","G21","G22","G23","G24","G25","G26","G29","G31","G32","G36","G44","G45","G50","G52","G55","G56","G59","G60","G61","G62","G66","G77","G79","G79","G80","G81","G85","G87","G89","G90","G91","G92","G93","G94","L01","L04","L05","L12","L18","L23","L25","L29","L32","L36","L44","L52","L59","L60","L61","L62","L66","L67","L79","L80","L81","L87","L89","L91","L92","L93","L94","M01","M04","M05","M12","M18","M21","M23","M24","M25","M29","M30","M32","M36","M44","M50","M52","M55","M56","M59","M60","M61","M62","M66","M67","M68","M79","M80","M81","M87","M89","M91","M92","M93","M94","Q12","Q18","S67","S74","S89")

list_rev <- c("B01","B21","B22","B30","B42","B46","B50","B59","B79","B80","B89","B91","B92","B93","B94","C21","C30","C42","C46","C50","C79","C80","C89","C91","C92","C93","C94","D21","D30","D42","D46","D50","D79","D80","D89","D91","D92","D93","D94","T01","T09","T10","T11","T12","T13","T14","T15","T16","T19","T20","T21","T22","T23","T24","T25","T27","T28","T29","T40","T41","T50","T51","T53","T99","A01","A03","A09","A10","A12","A16","A18","A21","A36","A44","A45","A50","A56","A59","A60","A61","A80","A81","A87","A89","U01","U11","U20","U21","U30","U40","U41","U50","U95","U99","A90","A91","A92","A93","A94","X01","X02","X05","X08","Y01","Y02","Y04","Y11","Y12","Y51","Y52")

# For empty entry replace by zero
for(i in 1:length(list_exp)){
  
  if(length(which(colnames(slgf) == list_exp[i]))==0){ 
    slgf$new <- 0
    colnames(slgf)[ncol(slgf)] <- list_exp[i]
  }  
}


for(i in 1:length(list_rev)){
  
  if(length(which(colnames(slgf) == list_rev[i]))==0){ 
    slgf$new <- 0
    colnames(slgf)[ncol(slgf)] <- list_rev[i]
  }  
}

#######
## Define variables
#######
slgf$total_taxes <- slgf$T01+slgf$T09+slgf$T10+slgf$T11+slgf$T12+slgf$T13+slgf$T14+slgf$T15+slgf$T16+slgf$T19+slgf$T20+slgf$T21+slgf$T22+slgf$T23+slgf$T24+slgf$T25+slgf$T27+slgf$T28+slgf$T29+slgf$T40+slgf$T41+slgf$T50+slgf$T51+slgf$T53+slgf$T99 

slgf$property_tax <- slgf$T01
slgf$sales_tax <- slgf$T09+slgf$T10+slgf$T11+slgf$T12+slgf$T13+slgf$T14+slgf$T15+slgf$T16+slgf$T19
slgf$income_tax <- slgf$T40
slgf$corporate_income_tax <- slgf$T41
slgf$vehicle_license <- slgf$T24
slgf$Other_taxes <- slgf$total_taxes - slgf$property_tax - slgf$sales_tax - slgf$income_tax - slgf$corporate_income_tax - slgf$vehicle_license

slgf$current_charges <- slgf$A01+slgf$A03+slgf$A09+slgf$A10+slgf$A12+slgf$A16+slgf$A18+slgf$A21+slgf$A36+slgf$A44+slgf$A45+slgf$A50+slgf$A56+slgf$A59+slgf$A60+slgf$A61+slgf$A80+slgf$A81+slgf$A87+slgf$A89

slgf$miscelaneous_gen_revenue <- slgf$U01+slgf$U11+slgf$U20+slgf$U21+slgf$U30+slgf$U40+slgf$U41+slgf$U50+slgf$U95+slgf$U99

slgf$utility_revenue <- slgf$A91+slgf$A92+slgf$A93+slgf$A94

###
# Aggregate by cities and school districts for graphing
###

slgf <- slgf[,c("type","utility_revenue","current_charges","miscelaneous_gen_revenue","property_tax","sales_tax","Other_taxes")]
slgf <- slgf[slgf$type %in% c(4,5),]
slgf$type <- ifelse(slgf$type ==4,"Cities","School districts")
slgf <- aggregate(slgf[,2:7],by = list(type = slgf$type),FUN = sum)

# Graphs are made in Excel from the file 
write.csv(slgf,"Output/fig A7.share_revenues.csv",row.names = FALSE)

rm(list = setdiff(ls(),c("df","city")))


###########################################
# Figure B1: Comparing the implied cost of bond vs. parcel levy referendums

###########################################


dfgraph <- df[df$type %in% c("Bond","Parcel Levy"),]
dfgraph$cost_comparison <- ifelse(!is.na(dfgraph$cost_bond_ph),dfgraph$cost_bond_ph,dfgraph$levy_parcel) 

q <- ggplot(dfgraph,aes(x = cost_comparison,group = type,color = type))+
  geom_density(size =0.8)+
  scale_x_continuous(name = "Imputed cost per parcel or housing units",label = dollar,limits = c(0,1000))+
  scale_color_manual(values=c("#B31B1B","#222222"),name = "Type of referendum:") +  
  theme_minimal()


ggsave(plot = q, file="Output/fig B1.bond_cost_vs_levy_density.png" ,bg = "white",width = 8,height = 5)
rm(list = setdiff(ls(),c("df","city")))

###########################################
# Figure C1: Income of school districts proposing referendums

###########################################

source("Rscripts/11bis.construct_matrix_of_refdums.R")

acs <- read.csv("Data/Raw data/ACS_SchoolDistricts_2010-2021.csv",header = TRUE)
acs <- cbind(LEAID = sprintf("%07d",acs$GEOID),acs)
acs$GEOID <- NULL
acs <- acs[acs$year ==2017,c("LEAID","income_mean")]

rdm <- merge(rdm,acs,by = c("LEAID"))

####
rdm$income_mean <- as.numeric(rdm$income_mean)

b <- aggregate(list(income = rdm$income_mean),by = list(refdum = rdm$refdum_binary,post = rdm$post),FUN = mean)
c <- aggregate(list(sd = rdm$income_mean),by = list(refdum = rdm$refdum_binary,post = rdm$post),FUN = sd)
d <- aggregate(list(nb = rdm$income_mean),by = list(refdum = rdm$refdum_binary,post = rdm$post),FUN = length)
b <- merge(b,c,by = c("refdum","post"))
b <- merge(b,d,by = c("refdum","post"))
b$refdum <- ifelse(b$refdum==1,"With Referendum","Without Referendum")

b$stdev <- (b$sd * 1.96) /sqrt(b$nb)

b$post <- ifelse(b$post==0,"Pre-TCJA","Post-TCJA")


q <- ggplot(data = b, aes(x = post, y = income, fill = post)) +
  geom_bar(stat = "identity", position = position_dodge(), alpha = 0.75)+
  facet_grid(. ~refdum)  +
  theme_minimal()+
  geom_errorbar(aes(ymin=income-stdev, ymax=income+stdev), width=.2,
                position=position_dodge(.9)) +
  scale_fill_manual(values=c("#B31B1B","#222222"),name = "School districts with:") +  
  scale_y_continuous(name = "Mean household Income",labels = dollar)+
  labs(x = "")+
  theme(legend.position = "none")

ggsave(plot = q, file="Output/fig C1.Mean_income.png" ,bg = "white",width = 8,height = 5)
rm(list = setdiff(ls(),c("df","city")))


